function priceapprox=price
global x r T t rho K % price parameters
global az bz mz % Zt, OU process parameters
global ay by my % Yt, OU process parameters
global eps v% = 1/ay
global sigma1sqhd sigma1sqbar
global Npath betaK

x = 0; r = 0; T = 0.5; t = 0; rho = 0.5; K=1; % x has to be nonzero in order P1 != 0 
az = 1;   bz = 1; mz=1; % let Zt be the simplest OU proess
ay = 1000; by = 1; my=1; % Yt 's coeff as OU process, ay = 1/eps;
eps = 1/ay; v = sqrt(by^2/2/ay); % page 68
sigma1sqhd = @(s) s.^2; % sigma1 = s
sigma1sqbar = barsigmasq; % <sigma1sq> = <s^2> = v^2+my^2 see page 68 page 86

Npath = 50000; betaK=20; 

eta=randn(Npath,betaK);
% etasquare=1/sqrt(2)*(eta.^2-1);
% etaeta=zeros(Npath, betaK*(betaK-1)/2);
% counter=1;
% for ii=1:betaK-1
%     for jj=ii+1:betaK
%         etaeta(:,counter)=eta(:,ii).*eta(:,jj);
%         counter=counter+1;
%     end
% end

intZtsq = getintZtsq(eta);
intZt = getintZt(eta);


% mu = x+r*(T-t)-Bk(j,T)^2/2*intZsq;
% sigma = Bk(j,T)^2*intZsq;

sampleprice0=zeros(Npath,1); sampleprice1=zeros(Npath,1);
[v2,v3] = coeffv3;

for ii=1:Npath
    [P0,ptptP0,ptptptP0] = conditionedprice(intZtsq(ii));
    sampleprice0(ii) = P0;
    v2bar = -v2*intZt(ii); v3bar = -v3*intZt(ii);
    sampleprice1(ii) = v2bar*x^2*ptptP0+v3bar*x^3*ptptptP0;
    disp(ii)
end

price0 = sum(sampleprice0)/Npath;
price1 = sum(sampleprice1)/Npath;

priceapprox  = price0+price1;
end

function Zt = getZt(u,eta)
global betaK T

s=analytichandle;

Zt = s.meanzt(u)+s.ztxi1(u,T)*eta(:,1);

for ii=2:betaK
    Zt = Zt + s.ztxik(ii,u,T)*eta(:,ii);
end

end


function intZt = getintZt(eta)
global betaK t T

s = analytichandle;
intZt = s.mean0(t,T)+s.xi1(t,T)*eta(:,1);
for ii=2:betaK
    intZt   = intZt+s.xik(ii,t,T)*eta(:,ii);
end
end

function intZtsq = getintZtsq(eta)
global t T Npath
time = linspace(t,T,10001);
intZtsq=zeros(Npath,1);
for ii=1:10001
    intZtsq = intZtsq + getZt(time(ii),eta).^2;
end
dt=(T-t)/100000;
intZtsq = (intZtsq-getZt(time(1),eta).^2/2-getZt(time(10001),eta).^2/2)*dt;


% s = analytichandle;
% intZsq = s.sqmean0(t,T)+s.sqxi1(t,T)*eta(:,1)+...
%     +s.sq1(t,T)*etasquare(:,1);            % construct int_t^T Z(u)^2du
% 
% intZ   = s.mean0(t,T)+s.xi1(t,T)*eta(:,1);
% 
% 
% for ii = 2:betaK             % project on Xi
%     intZsq = intZsq+s.sqonxik(ii,t,T)*eta(:,ii);
%     intZ   = intZ+s.xik(ii,t,T)*eta(:,ii);
% end
% 
% for ii=2:betaK               % project on xi square
%     intZsq=intZsq+s.sqonxisq(ii,t,T)*etasquare(:,ii);
% end
% 
% for ii=2:betaK               % project on 1_xi
%     intZsq=intZsq+s.sqmix1k(ii,t,T);
% end
% 
% counter = 1;
% for ii=2:betaK-1             % project on xi_i xi_j
%     for jj=ii+1:betaK
%         intZsq=intZsq+s.sqmixij(ii,jj,t,T)*etaeta(:,betaK-1+counter);
%         counter=counter+1;
%     end
% end

end


function [P0,ptptP0,ptptptP0] = conditionedprice(intZsq)
global x r T t K sigma1sqbar

mu = x+r*(T-t)-sigma1sqbar/2*intZsq;
sigma = sigma1sqbar*intZsq;

P0 = quadgk(@(s) max(s-K,0).*lognpdf(s,mu,sigma),-Inf, Inf);
ptptP0 = quadgk(@(s) ((mu-log(s)).^2/sigma^4-1/sigma^2).*max(s-K,0).*lognpdf(s,mu,sigma),-Inf, Inf);
ptptptP0 = quadgk(@(s) ((mu-log(s)).^3/sigma^6+3*(mu-log(s))/sigma^4).*max(s-K,0).*lognpdf(s,mu,sigma),-Inf, Inf);

P0 = exp(-r*(T-t))*P0;
ptptP0 = exp(-r*(T-t))*ptptP0;
ptptptP0 = exp(-r*(T-t))*ptptptP0;
end

function [v2,v3] = coeffv3
global rho ay v

phiprime = @newapproach;
average=quadgk(@(y) phiprime(y).*y, -Inf, Inf); 
% here phi' is actually phi'(y)*normpdf(y,m,v), so average actually
% = y*phi'(y)*normpdf(y,m,v)=<yphi'(y)>

v3 = v/sqrt(2*ay)*rho*average;
v2 = 2*v3;
end

function result=newapproach(y) % the verification see ../testing/testcoeff3.m
global v my sigma1sqbar
z=y-my;
%A1=@(z) (z<0).*(sqrt(pi)/4*erfc(-z)-z/2.*exp(-z.^2))+(z>=0).*(sqrt(pi)/2-z.*exp(-z.^2)/2-sqrt(pi)/4*erfc(z));
A1=@(z) sqrt(pi)/2-z.*exp(-z.^2)/2-sqrt(pi)/4*erfc(z);
result = 1/sqrt(2*pi)/v*(sqrt(2*v^2)^3*A1(z/sqrt(2)/v)-2*v^2*my*exp(-z.^2/2/v^2)+my^2*sqrt(2*pi*v^2)/2*erfc(-z/sqrt(2*v^2)))-sigma1sqbar/2*erfc(-z/sqrt(2*v^2));
result = result/v^2;
% result <==> quadgk(@(s) (sigma1(s)-sigma1sqbar).*normpdf(s,my,v), -Inf, y)      /v^2;
% which is phi'(y)*normpdf(y,m,v);

end

function y = barsigmasq
global my v %sigma1sqhd
%y = quadgk(@(s) sigma1sqhd(s).*normpdf(s,my,v), -Inf, Inf);
y = v^2+my^2;
end